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ABSTRACT 


The threat of collision between an asteroid or a comet and the Earth has been well 
documented. Mitigation of such a threat can be accomplished by destruction of the threat 
or by perturbing the threat object into a safe orbit. Following a summary of proposed 
mitigation techniques, this thesis investigates the impulse required to safely perturb a 
threatening Earth-Crossing Asteroid (ECA). While previously published analysis 
included only two-body approximations to the impact geometry, this thesis adds the 
effect of the Earth's gravitational field to more closely approximate reality. The results 
indicate that third-body effects are strongest on ECA's in a nearly circular heliocentric 
orbit, where the minimum required AV can be several times larger than that calculated 
using two-body approximations. To determine the minimum AV required for mitigation, 
MATLAB? 's sequential quadratic programming (SQP) algorithm is applied to a 
constrained optimization problem. Third-body effects were added to a previously 
published two-body optimization by modifying the boundary conditions. With 
knowledge of the minimum AV requirements, the capability of current impulsive 
mitigation technology is analyzed. For asteroids of median density in co-planar orbits, a 
single 24 Mt nuclear explosive impulse applied earlier than 3 years before impact can 
effectively mitigate a threat with a diameter of 6 km. The capability significantly 


decreases with shorter warning times. 
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L INTRODUCTION 


Who knows whether, when a comet shall approach this globe to destroy it, 
as it often has been and will be destroyed, men will not tear rocks from 
their foundations by means of steam, and hurl mountains, as the giants are 
said to have done, against the flaming mass? - And then we shall have 
traditions of Titans again, and of wars with Heaven. 

- Lord Byron, 1822 (Medwin, 1824, p. 185) 


Although steam is no longer the motive force of choice, the defense that Lord 
Byron dreamed of 150 years ago is finally possible today. For the first time in the 
history of humanity, the technology exists which could potentially protect the Earth 
against an impact from an asteroid or comet. In addition, the threat of impact from a 
cosmic body that Lord Byron addressed has, within the last 20 years, become widely 
accepted. However, while the defensive technology exists, an active planetary defense 
plan does not. At this incipient stage in researching impact hazards, there are still many 
challenges to be met. 

Before considering the technical challenges to a defense plan, one must first 
answer the question: what is the threat? Cosmic bodies in orbits that bring the 
possibility of impact with the Earth are commonly termed Near Earth Objects (NEOs). 
Generally, NEOs are comets or asteroids that have a perihelion radius of less than 1.3 
AU. These are orbits that have the potential to be perturbed by the gravitational force of 
the Earth and other planets into orbits that impact the Earth (Cheng, 1994, p. 651). | 

Lord Byron's quotation refers to comets, which are distinguished from asteroids 
by outgassing activity or the existence of a tail (Cheng, 1994, p. 651). Comets are 
believed to be composed of ice and rock (Elder, 1997, p. 11), and can approach in such 
highly eccentric orbits that they can impact the Earth at velocities greater than 50 km/s 
(Morrison, 1994, p. 61). The composition of asteroids varies widely from metal to the 


same ice and rock found in comets (Morrison, 1994, p. 61). Because of the lower 


eccentricity of their orbits, asteroids generally impact the Earth at velocities on the order 
of 20 km/s (Morrison, 1994, p. 61). 

The technical challenges to a NEO defense plan can be grouped into two 
categories: hazard detection and hazard mitigation. 

Current theory suggests that only a small percentage of all potentially hazardous 
Near Earth Objects (NEOs) have been discovered. It is believed that all Earth-Crossing 
NEOs with diameters of 6-12 km have been discovered. With smaller diameters, the 
detection completeness decreases dramatically: about 35% with diameters of 3-6 km, 15% 
with diameters from 2-4 km, and 7 % with diameters from 1-2 km. This low detection 
completeness implies that a significant number of objects with the potential of causing a 
global disaster have not yet been detected. Some specific challenges for detection are: 
correcting a bias in discoveries that results from the majority of searches being conducted 
in the Earth's northern hemisphere, improving discovery follow-up to achieve a reliable 
orbit calculation (Carusi, 1994, p. 145), decreasing the minimum possible detection 
magnitude, and coherently processing detection data (Bowell, 1994, p. 194). 

Once detected, the orbital parameters of the NEO must be accurately determined. 
Errors in the estimation of a NEO's orbital parameters can propagate over time to produce 
a drastically inaccurate estimate of impact probability or time. A thorough understanding 
of this error propagation is critical to the planning of a defense mission. Another 
challenge is to model the non-gravitational effects in comets caused by rocket-like 
outgassing (Yoemens, 1994, p. 257) so that accurate impact probabilities can be 
calculated. 

This thesis concentrates on some of the challenges to the mitigation of a NEO 
hazard. Chapter II presents a detailed study of mitigation ideas. The thesis then 
continues with the goal of answering two questions: 


1. What is the minimum AV required to deflect a threat-NEO into a safe orbit? 


2. What is the maximum NEO size that can be deflected with a single defense 


mission using current technology? 


The challenges addressed in this thesis are only a fraction of those that must be 


overcome before a Planetary Defense system can be fielded. 





II. STATEMENT OF THE PROBLEM 


A. HAZARD MITIGATION TECHNIQUES 


All proposed NEO hazard mitigation techniques (Simonenko, 1994; Wood, 1994; 
Ahrens, 1994; Solem, 1994; Canavan, 1994; Meissinger, 1995; Melosh, 1994; Shafer, 
1994) attempt to accomplish one of two goals: either destroy the NEO, producing 
fragments that will not inflict damage if they impact the Earth, or deflect the NEO into an 
orbit that is no longer hazardous. Both of these techniques require a controlled coupling 
of energy to the NEO. Much of the analysis of hazard mitigation techniques thus far has 
revolved around determining the best source of energy for hazard mitigation (Canavan, 
1994; Ahrens, 1994; Simonenko, 1994; Shafer, 1994; Solem, 1994; Willoughby, 1994; 
Melosh, 1994). 

As with all space missions, the mass of the payload is a critical driver in the 
selection of energy source. Table 2.1 gives a comparison of the specific energies of many 
potential energy sources. If energy source selection were based solely on payload mass, 


nuclear explosive would obviously be the best choice. 


High Explosive 6 MJ/kg 
Kinetic Energy (10 km/s) 50 MJ/kg 


Nuclear Explosive 4*10% MJ/kg 





Table 2.1 Specific Energies (after Shafer, 1994) 


It is instructive at this point to discuss terminology used in describing nuclear 
explosives. In most technical applications, the size of a nuclear explosive is expressed as 
the amount of energy the explosive can produce. This energy is usually measured in 
terms of an equivalent mass of TNT that would produce the same explosive power. 
Using this convention, then, 1 MT = 4.2 x 10%” Joules (Morrison, 1994, p. 61). Table 2.2 


presents the mass of a nuclear explosive charge based on its yield. 


1 Mt 0.5 ton 


10 Mt 3 to 4 ton 
100 Mt 20 to 25 ton 





Table 2.2 Yield Versus Mass for Nuclear Explosive Devices (after Simonenko, 1994, p. 
931) 


However, mass is only one of several factors that must be considered in choosing 
an energy source. The following is a compilation of proposed NEO hazard mitigation 


techniques, based on the mitigation goal. 


1. Fragmentation and Dispersal 


If sufficient energy is delivered to the NEO to fragment it, the hazard could be 
mitigated in several ways. First, differential velocities in the individual fragments of the 
destroyed NEO would result in a dispersed debris cloud. Some small scale fragmentation 
experiments suggest that if a NEO were destroyed 1 orbit prior to impact, as little as 1% 
of the original NEO mass would impact the Earth (Ahrens, 1994, p. 919). 

The second mitigation effect of fragmenting the NEO is to reduce the size of the 
meteor that reaches the Earth. In order to be effective, a mitigation effort would be 
required to reduce the impactor size to the point that it would no longer generate damage 
on impact with the Earth. The atmosphere can provide a great deal of protection from 
some meteors, but its effects are highly dependent on the impact velocity and the 


composition of the NEO. Figure 2.1 shows the relation between the composition of a 


NEO, its impact velocity, its size, and the amount of destruction that can result following 


atmospheric ablation and impact. 
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Figure 2.1 The radius of destruction around the impact point due to the atmospheric blast 
wave (from Hills, 1993, p. 1133). 


But what is the best way to fragment the NEO? For the mass reasons discussed 
above, nuclear explosives are a very attractive option for the fragmentation energy source. 


Smaller, more loosely bound NEOs could be destroyed by a surface explosive. However, 


larger and more solid NEOs would require that the explosive be buried. Burial would 
provide significantly greater coupling of energy to the NEO (Ahrens, 1994, p. 917). 

The requirement for charge burial adds some complexity to the mitigation effort. 
Wood (1994) suggests some ways to effect this burial. A penetration depth of 
approximately 10 meters could by accomplished by use of hyper-velocity (kinetic 
energy) blasting prior to arrival of the nuclear charge. For larger NEOs, a string of 20-100 
0.5 MT nuclear explosives could be used for penetration. Once the desired depth is 
achieved, a 20-300 MT charge can be directed into the void and detonated to achieve final 
NEO destruction. 


Ahrens and Harris (1994, p. 922) estimate that a charge buried at optimal depth 


can destroy a hazardous NEO as summarized in Table 2.3: 


NEO Diameter Required Nuclear Explosive Yield 





Table 2.3 Nuclear Explosive Yield required to fragment hazardous NEOs (after Ahrens, 
1994, p. 922) 


While analysis indicates that nuclear explosives can successfully fragment a NEO, 
there are hazards associated with this energy source. First 1s the concern, presented by 
Meissinger (1995, p. 16), that the NEO fragments which reach the Earth following nuclear 
destruction could be radioactive. The most significant problem with using nuclear 
weapons in NEO mitigation is the tremendous political, sociological, economic, and 
management aspects of the weapons. 

The alternatives for NEO fragmentation use the kinetic energy inherent in the 
closing speeds between the NEO and an interceptor. Wood (1994) proposes a "Jack- 
Hammer" rubblization technique in which waves of hyper-velocity projectiles are directed 


towards the threat object. Each projectile vaporizes a small portion of the object, leaving 


a hole for the next projectile to enter. The cumulative effect of the impact of many 
projectiles is several deep holes and a somewhat fragmented NEO. 

Woods’ team further refines this concept (Teller, 1995) into Hypervelocity 
Projectile Array-Sheets and Hypervelocity Projectile Lattices. The Sheets are a two 
dimensional plane with multiple projectiles connected to each other. The connection 
helps to ensure that follow-on projectiles will be directed to impact spots of previous 
projectiles. The Lattice combines several Sheets into a three dimensional pulverization 
effort. | 

Thus, there are several options for fragmentation and dispersal of a threat-NEO. 
There are, however, arguments against using fragmentation as a mitigation goal. In order 
for fragmentation to significantly reduce the amount of NEO mass that eventually 
impacts the Earth, the NEO must be destroyed early. Thus, Ahrens and Harris conclude 
that "...fragmentation is likely to be a safe choice only for long lead-time response 
(decades) or for relatively small bodies where the fragments may be allowed to hit the 
Earth (Ahrens, 1994, p. 919)." If the target were pulverized late, it may make the 
problem worse. The resulting fragments would present multiple targets to any follow-on 


mitigation effort (Solem, 1994, p. 1028). 


D Orbital Deflection 


Many NEO orbital deflection schemes work by accelerating mass away from the 
NEO, thus changing its momentum. Chemical rockets, which have seen extensive use in 
space for changing the momentum of man-made satellites, are only appropriate for NEOs 
with diameters of less than 100 meters (Wood, 1994, p. 11). The mass of fuel required 
to rendezvous (by matching velocity) with the threat and then to displace it would be 
otherwise be excessive. 

Another option is to use a mass driver. A mass driver uses electromagnetic forces 


to accelerate buckets containing material from the surface of the NEO (Melosh, 1994, p. 


1117). Thus, a mass driver has an advantage over chemical rockets in that it avoids the 
need to transfer propellant mass from the Earth. While the technology for mass drivers 
has been available since the mid-1970's (Melosh, 1994, p. 1117), they are nonetheless 
very complex. Meissinger (1995) points out that mass drivers would require large-scale, 
complex robotic operations for construction and processing of the NEO's soil. Also, the 
mass driver concept is critically dependent on the soil conditions of the NEO (Meissinger 
1995). At this point, little is known about the composition or soils characteristics of 
asteroids and comets. 

Just as the kinetic energy of the threat-NEO was used in NEO destruction 
schemes, there are schemes to use this kinetic energy for NEO deflection. Interestingly, 


for kinetic energy mitigation schemes, Canavan (Canavan, 1994, p. 102) suggests that 


. m : 
",..the mass required scales on —, so faster NEOs present less of a threat because of their 
V 


higher specific energy." In these schemes, a projectile impacts the threat-NEO and the 
crater ejecta combines with the momentum of the projectile to produce an impulse. 
Analysis conducted by Ahrens and Harris (1994, p. 904) suggests that the energy from a 
kinetic impact is coupled much more efficiently than the energy from nuclear weapons 
detonated at the surface of the NEO. Figure 2.2 shows an estimate of the capability of 
kinetic energy deflectors. The three lines show the capability of an impactor with a 
diameter of 1m, 10m, and 100m. Melosh (Melosh, 1994, p. 1115) even suggests a 
"billiards shot" scenario, in which a small asteroid is displaced into a larger asteroid. 
Wood (1994) also refines the kinetic energy approach into a "hypervelocity sand-blaster." 
In this scheme, a steady stream of projectiles are directed at the threat-NEO, combining 


their effect to produce the required momentum change. 
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Figure 2.2 Capability of Kinetic Energy Deflectors (from Melosh, 1994, p. 1116) 


The next potential technique for orbital deflection is the use of nuclear explosives. 
Momentum change can be imparted on a threat-NEO using a nuclear explosive in one of 
two fundamental ways: a stand-off radiative explosion or a surface explosion. For the 
stand-off case, the explosive is DE to have a substantial fraction of its yield as 
energetic neutrons and gamma rays. These rays irradiate the surface of the NEO that is 
exposed to the explosion. The irradiated surface then expands and ablates away from the 
NEO. As the material ablates away, the momentum of the NEO changes. The surface 
explosion works in much the same way as a kinetic defense scheme by creating a crater. 
The ejecta departing from this crater creates a momentum change in the NEO. (Ahrens, 
1994, p. 910). 

Not surprisingly, it is not yet clear which nuclear explosive method has the 
greatest mitigation capacity. Surface nuclear explosions are generally considered more 
effective in coupling energy into the NEO. Canavan (Canavan, 1994, p. 105) concludes 
that "...the energy required for stand-off is greater than that for slightly subsurface bursts 
by a factor of about 40." Solem (1994, p. 1032) reaches the same conclusion, stating that 
a stand-off deflection would require a substantial increase in interceptor mass over a 


surface burst. Ahrens and Harris (1994, p. 917), in contrast, propose that "...surface 


11 


explosions appear to be not substantially better than radiative stand-off explosions, in 
deflecting NEOs." One significant disadvantage to the surface explosion is the risk of 
inadvertently fragmenting the threat object. As pointed out previously, a fragmented 
NEO presents multiple targets for any follow-on mitigation effort. But there is fairly 
consistent agreement that the interceptor weight for a nuclear mitigation effort is several 
orders of magnitude less than that of a kinetic energy effort (Solem, 1994, p. 1032). 

Lasers also have a potential use in NEO hazard defense. With this technique, a 
high energy laser is directed at the surface of the threat-NEO. This creates a thermal flux 
which ablates the surface, much like the radiative heating does in stand-off nuclear 
explosions. One estimate suggests that a laser output of 1 GJ/s for 12 uninterrupted days 
could match the energy fluence of a 1 Mt nuclear burst. (Shafer, 1994, p. 965) 

Solar sails have also been proposed. Solar sails use solar radiation pressure to 
provide a motive force. This force is small, but consistent over long periods of time. The 
long build up could provide a significant change in momentum of a hazardous NEO. 
However, Melosh points out that "...truly enormous structures are necessary to deflect 
asteroids in the 1 to 10 km diameter range." Figure 2.3 demonstrates the required size. 
The lines are for solar sails of 10 km, 100 km, and 1000 km diameters. The technology 


does not yet exist to build such enormous solar sails at great distances from the Earth. 
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Figure 2.3 Capability of Solar Sails (from Melosh, 1994, p. 1120) 


One of the more interesting non-nuclear approaches is a solar collector. Melosh 
(Melosh, 1994, p. 1120) has done some extensive analysis of this approach and concludes 
that, for non-nuclear alternatives, it is "...an approach that is arguably better than any 
other previously proposed.” This scheme also uses a solar sail, but instead of using the 
solar sail directly as a motive force the sail is used to focus sunlight onto the surface of 
the NEO. The surface 1s thus vaporized, and the ablation of surface material produces a 
momentum change. Melosh indicates that a 1 km diameter solar collector, which can be 
launched by the space shuttle, can deflect asteroids up to 3.4 km in diameter after 


operating for a year (Melosh, 1994, p. 1125). 


B. PREVIOUS ORBITAL DEFLECTION ANALYSIS 


The fundamental goal of this thesis is to investigate the effects of third-body 
gravitation on orbital deflection requirements. Previous research has simplified the 
problem by assuming two-body orbital mechanics between the Sun and the threat-NEO. 
This assumption neglects perturbations due to the Earth's gravity. While these 
perturbations may not be present until the terminal phase on the impact scenario, they 


affect both long and short warning time analyses as described below. 
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1. Long warning times 


Most analyses of long warning time threats have concentrated on changing the 
phasing of the NEO. Ahrens (1994, p. 903) approximated the required AV for deflection 
by comparing the position of the perturbed NEO orbit with that of the unperturbed orbit. 
He then calculates the velocity change required to produce a displacement of 1 Earth 
radius (R,). Investigation of third body effects and patched conic theory indicates that a 
displacement of 1 R, may be insufficient. 

Earlier work done at the Naval Postgraduate School (Knudson, 1995; Park, 1997; 
Elder, 1997) refines the estimates of the required AV. All of these analyses, however, 
begin with a two-body approximation. One would like to know how much third-body 


gravitation will affect these results. 


2: Short warning times 


Ahrens (1994, p. 901), Meissinger (1995), and Solem (1994, p. 1015) all assume 
rectilinear motion for short warning time calculations. Once again, this approximation 
neglects the Earth's gravitational effects, which will cause the NEO to deflect towards the 


Earth in a hyperbolic orbit. 

A brief look at hyperbolic orbits will give an idea of how Earth's gravity will affect 
rectilinear approximations. When a sideways motion is imparted such that a 
(gravitationally free) miss distance of 1 R4 is achieved, the result is to establish a 
hyperbolic orbit with a semiminor axis (b) of 1 Rg. Using a typical impact velocity of 20 


km/s, actual perigee radius (equation from Brown, 1992, p. 27) can be calculated as 


follows: 
2 
H u 2 
r = -— + | — | +b = 0.856R 
Ge SI 8 
where 
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Və = impact velocity 
u = Gravitational Constant of the Earth 
Thus, preliminary approximations of energy required to deflect a threat may result 
in errors of 14% in miss distance, in the terminal case. 
This thesis intends to further investigate the Earth's gravitational effects on energy 


required to defend against an impacting NEO. 
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lil. PROBLEM FORMULATION 


A. PROBLEM STATEMENT 


Given a NEO with an orbit that confirms an impending collision with the Earth, 


and given that orbital deflection is the chosen mitigation goal, it is necessary to know the 


AV (both magnitude and direction) required to deflect the asteroid into a safe orbit. 


Determination of this required AV must include multi-body gravitational effects. A 


detailed understanding of third-body perturbations will help to refine previous two-body 


approximations. 


The goal of this thesis is to investigate the effects of the Earth's gravitation on the 


AV required to deflect an Earth-impacting NEO. 


B. ASSUMPTIONS 


Several assumptions were made to simplify the analysis. 


The Earth and the threat-NEO are in co-planar orbits. 

The Earth is considered to be in a perfectly circular heliocentric orbit at a radius 
of 1 AU. 

The threat-NEO is originally in an elliptical heliocentric orbit. NEOs in 
hyperbolic or parabolic heliocentric orbits were not considered. 

The control maneuvers are impulsive. 

Other than the impulse, no non-gravitational forces (such as outflowing of gas 


and dust in comets) are included. 


C. HYPERBOLIC MAPPING AND ANALYSIS 


The first step was to understand how a heliocentric elliptical orbit maps into a 


geocentric hyperbolic orbit. The primary goal was to determine the velocity of the NEO 
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with respect to the Earth (Vreo ) at Earth impact. This required some vector analysis 


combined with orbital mechanics and geometry. 


Once an equation for Vreo was found, an understanding of the dynamics can be 
e 


aided by finding its maximum. The maximum was determined in two ways. First, 


Mathworks' MATLAB*^ was used to numerically create a 3-D plot of Veo, versus the 


NEO's heliocentric semimajor axis (a ) and its heliocentric eccentricity (e ). To confirm 
the numerical result, Warterloo's symbolic manipulator, Maple V™, was used. The 


original equation for Vreo, was maximized analytically. The plot was identical to that of 


the MATLAB* analysis, thus verifying both analyses. 


D. PATCHED CONIC APPROXIMATION 


Once the hyperbolic mapping was realized, the analysis was extended for use in a 
patched conic approximation. The advantage of using this approximation was that 
previous two-body analysis could be adapted to account for third body effects. 

The most useful aspect of the patched conic approximation, for this analysis, was 
the B-plane and the impact radius (b,). Figure 3.1 illustrates these parameters. The B- 
plane is orthogonal to the asymptote of the approach hyperbola and placed at a large 


distance from the Earth. The impact radius (b,) is the semiminor axis of a hyperbola that 


has a perigee equal to the radius of the surface of the Earth. 


Trajectory 
Plane 


impact Radius, bi 
Target Point, b 





Figure 3.1 B-plane and Impact Radius (from Brown, 1992, p. 117) 


The impact radius is a function of Vreo , which was solved for in the hyperbolic 


analysis. Once this radius is determined, it can replace previous miss distances of 1 R, in 


two-body analysis. Thus, third-body effects are accounted for while maintaining the 


relative simplicity of two-body analysis. 


E. NUMERICAL TRAJECTORY OPTIMIZATION 


Park, Elder, and Ross (Park, 1997) use MATLAB? 's sequential quadratic 
programming method to numerically solve a constrained optimization problem to 
determine the minimum AV for deflecting an Earth-crossing asteroid. As referred to 
earlier, this code used a two-body approximation. It is possible to include the effects of 
third-body gravitation in this code by adding the impact parameter in two different ways. 

Discussing these modifications first requires a brief description of the code. The 


code requires a target separation distance. This target separation 1s the minimum 
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allowable separation between the Earth and the perturbed NEO. Following an impulse, 
the Earth's orbit and the NEO's orbit are propagated using the two-body equations of 
motion to determine the minimum separation. The optimization problem minimizes the 
impulse by equalizing the propagated and the target separations. 

The first modification of the code is to increase the target separation distance. If 
this target separation is increased to a minimum impact radius in the B-plane, then the 
Earth's gravity will not cause an impact. 

Second, the constraints of the optimization problem can be modified. In its 
original form, the constraints require that the minimum separation equal some fixed 
distance. With the first modification, this fixed distance is the impact radius of the NEO's 
original orbit. Once the NEO's orbit changes (after application of the AV), there is a new 
impact radius. If this new impact radius is larger than the original impact radius, the NEO 
will still impact the Earth. The solution is to modify the constraint to include the impact 
radius directly. This modification requires that the minimum separation be equal to the 


impact radius. 
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IV. HYPERBOLIC ORBIT ANALYSIS 


A. CONIC SECTIONS 


Orbits in an inverse square gravitational field take the shape of conic sections: 
circular, elliptical, parabolic, or hyperbolic. Circular orbits are a special case of the 
elliptical orbit. These two types of orbits are defined when the object's energy is 
insufficient to escape the gravitational attraction of the central mass. The parabolic orbit 
is a transition between elliptical and hyperbolic orbits. An object passing a central mass 
on a parabolic trajectory would reach an infinite distance from the mass, but with zero 
velocity. Any object with a velocity less than that of a parabolic orbit will not escape the 
central mass, and thus will be in an elliptical orbit. Conversely, an object with a greater 
energy than that of a parabolic orbit will escape the central mass. When the object is 


traveling faster than this escape velocity, it is in a hyperbolic orbit. 


B. ELLIPTICAL ORBITS 


Before they become a threat to Earth, nearly all NEOs are established in an 


elliptical orbit about the Sun. Figure 4.1 shows the geometry of an elliptical orbit. 


Line of Apsides 


Apoapsis 





Figure 4.1 Elliptical Orbit (after Brown, 1992, p. 9) 


For the two-dimensional case considered in this thesis, only the semi-major axis 


(a), eccentricity (e), and true anomaly (v) are required to define an orbit and position. 
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a = semimajor axis, the distance from the center of the ellipse to the long edge 

e = eccentricity, defines the shape of the orbit (0 <e < 1 for an elliptical orbit) 

v — true anomaly, defines the position within the orbit 

The most important elliptical orbit relation required for this analysis is the NEO's 
heliocentric velocity as it reaches the orbit of the Earth. This velocity is given by 
(Brown, 1992, p. 19) 


2 
Y (4.1) 
r a 


Also useful is the distance of the NEO from the Sun as a function of the true 
anomaly, given by (Brown, 1992, p. 18) 
a(1 - e?) 


ee 4.2 
1+ecosu e, 


C. HYPERBOLIC ORBITS 


In an impact scenario, a NEO which was in an elliptical orbit about the Sun will 
transition, as it approaches the Earth, to a hyperbolic orbit about the Earth. An 
understanding of hyperbolic orbits, and this transition, is therefore fundamental to studies 
of the orbital mechanics of NEO impacts. 

The orbital parameters of a hyperbolic orbit, shown graphically in Figure 4.2, are 
similar to the familiar elliptical orbit parameters (Brown, 1992, p. 21): 

r, ^ periapsis radius 

a — absolute magnitude of the semimajor axis, the distance from the asymptote 

focus to the periapsis 

b — semiminor axis, distance from the asymptote to a parallel passing through the central boc 

c — the distance from the asymptote focus to the center of mass 


e = eccentricity = Y, (greater than 1 for hyperbolic orbit) 


ß = angle of the asymptote 
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6, = true anomaly of the asymptote 


ol 
A 


b c 


Line of Apsides 


\ 





Figure 4.2 Elements of a hyperbola (after Brown, 1992, p. 22) 


All of the orbital parameters of a hyperbolic orbit can be determined with a 
knowledge of two parameters: the velocity at an infinite distance from the central mass 
(V.,.) and the semiminor axis (b). From knowledge of V.,, the semimajor axis can be 


calculated by the relation (Brown, 1992, p. 28) 

V zye (4.3) 
a 

Once the semimajor axis (a) is determined, the eccentricity can be calculated by 


(Brown, 1992, p. 27) 
ee (4.4) 


Periapsis radius (r, ) is expressed as (Brown, 1992, p. 27) 
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Ask =; E +b? (4.5) 


And finally, the angle of the asymptote ( B ) is (Brown, 1992, p. 27) 


tan p — 2 (4.6) 
a 


D. PATCHED CONIC METHOD 


The patched conic method was originally developed to simplify the planning of 
interplanetary trajectories, yet it applies directly to analysis of NEO impacts on the 
Earth. 

Each celestial body, due to its mass, generates a gravitational field which affects all 
other celestial bodies according to Newton's Laws of motion and gravity (Wiesel, 1997, p. 
35) 

mr = y I (4.7) 
mi P L 

This is known as the N-body problem. An exact solution for more than two 
bodies does not exist (Weisel, 1997, p. 35). However, the patched conic method makes 
some minor approximations that reduce the N-body problem to solvable two-body 
problems. 

The patched conic method assumes that an object is influenced by the 
gravitational field of a planet only when it is within the planet's "sphere of influence." 
Beyond the sphere of influence, the object is considered to only be affected by the Sun's 
gravitation. The radius of the sphere of influence is somewhat nebulous, but Laplace 
suggests (Brown, 1992, p. 97), 

Mus Y^ 
Ro de) (4.8) 


sun 
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where 
Ryo, = radius of the sphere of influence of the planet 


R = mean orbital radius of the planet 


M 


planet 


= mass of the planet 


M 


Sun 


= mass of the Sun 


For Earth, this equates to a sphere of influence of 0.9 x 10°km. 
In the case of a NEO impact scenario, the NEO begins in an elliptical orbit about 
the Sun. Once within the sphere of influence of the Earth, the NEO's motion is described 


by two-body orbit equations for a hyperbolic orbit about the Earth. 


E. DETERMINING V, 


To define the orbit of the NEO as it approaches the Earth, it is necessary to 
determine V,. Figure 4.3 shows the vector geometry of an impact scenario, where: 


Vreg = the velocity of the NEO with respect to the Sun 

© 

Vo, = the velocity of the Earth with respect to the Sun 
© 


V wo 


Ye = the velocity of the NEO with respect to the Earth, defined as Va in this 


thesis 


Strictly speaking, Vreo, # V_, because the NEO is not at infinity with respect to 


the Earth. However, in the patched conic method the approximation is made as 


Vaso = Vo 
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Earth’s Orbit 
NEO’s Orbit 


X (Local Horizontal) 





Figure 4.3 Vector Geometry at Impact 
Using vector algebra, V4 can be expressed as: 


Voy. + Vyeo = ME (4.9) 
V = V = V, = V 4.10 
veo = Va = Voy 7 Vs, (4.10) 
The velocity of the Earth with respect to the Sun (Voy ) can be determined at any 
© 


time, using the equation for velocity in a circular orbit 


ao Pog (4.11) 
% 
with 
roy = the distance of the Earth from the Sun = 1 AU 
© 
To find the other unknown, Vago, , Some approximations are required. The 


© 
radius of the sphere of influence is very small in comparison to the radius of the Earth's 
orbit about the Sun. Thus, calculating the velocity of the NEO as it crosses the Earth's 


orbit is a very close approximation to the velocity as it crosses the Earth's sphere of 


26 


influence (SOI). In the worst case geometry, this approximation introduces an error in 


radius of 


0.9 x 10° 
"ioi x 100% = 222% 100% = 0.6% (4.12) 
PA 1.496 x 10° km 


Strictly speaking, this approximation changes equations 4.9 and 4.10 from 


equalities to approximations. 


Now the magnitude of MI is calculated by 
© 


2H Ve 
Y a NEO, 


V 


aA s 








(4.13) 


with 
a =the semimajor axis of the NEO with respect to the Sun 


For the direction of Vigo, , the flight path angle (y ) is required. The flight path 
© 


angle ıs a function of the position of the object within its orbit. This position is defined 
by the true anomaly ( v), shown in Figure 4.4. In an impact scenario wıth a known NEO 


orbit, the true anomaly at impact can be calculated as follows (where e : = the eccentricity 


of the NEO's heliocentric orbit): 


a(1- e?) 


Us = lay = (4.14) 
7% Jusen % mme VEO, 
© / IMPACT 


After manipulation, Equation 4.14 yields 


zelt 


( S (4.15) 
"we, » KEN | en, e| l 
Zo 


D 


Earth’s Orbit 


NEO’s Orbit 





Figure 4.4 Definition of true anomaly 


Given the position of the NEO and it's heliocentric eccentricity, the flight path 
angle is calculated by (Brown, 1992, p. 18) 
esinv 
N 


any = — (4.16) 
l + €COS Ung 
© 


Thus Equations 4.13, 4.15, and 4.16 combine to give both the magnitude and 


direction of Vio, . Returning to Equation 4.10, V_ can finally be determined: 
O 


= a 4.10 
V. MIT Vey ( ) 
V, = (Vaso cosy X + Varo siny? ) - Vy X (4.17) 

% © % 
\ Te, a Teg \ Tee a 
€ Sin U, go 
nup PC (4.18) 


1+ €C0S Ungog 


Digg ) = cos 
( 96 ) mapact | 
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Equations 4.18 are an important result. Although they look complicated, it is 
apparent that V, is a function of only two variables: the NEO's heliocentric eccentricity 
(e ) and semimajor axis (a ). The flight path angle and radius are set by the geometry of 
the impact. Since the Earth is assumed to be in a circular orbit, and the problem is planar, 


this geometry is defined by the NEO's orbit. 


F. MAXIMUM V, 


As previously discussed, V , is a critical value in defining the hyperbolic orbit of 
the NEO about the Earth's center of gravity. Now that V. has been directly related to e 
and a, a better understanding of this relation can be achieved by attempting to find its 


maximum. 


1. Variable Bounds 


If it is possible to find bounds for each parameter in equations 4.16, it may be 
possible to determine a maximum V.. Such a study will potentially provide an insight 
into the full geometry of the problem. 


To determine a value for the minimum semimajor axis of the NEO (amin), recall that 





a = — (4.19) 


A minimum value for the NEO's semimajor axis is mandated by the fact that the 
orbit must intersect the Earth's orbit for an impact to occur. This sets the value of r, to 
the radius of the Earth's orbit. The NEO also cannot have a perihelion radius less than the 
radius of the Sun (R,). Thus, amin is defined by 


p., 
AA 


(4.20) 
a 2 
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The maximum value of the semimajor axis is infinite. This is because the aphelion 
radius is unbounded, although the perihelion radius cannot be larger than the radius of the 


Earth's orbit. This becomes more clear when using the following elliptical relation 





a= (4.21) 


In an impact scenario with a NEO in an elliptical heliocentric orbit, the 
eccentricity varies as 0 <e <1. Thus, as the eccentricity approaches 1, the semimajor 


axis approaches infinity. The bounds of the semimajor axis can be summanized as 


RS 


The bounds of the NEO's eccentricity are set by the initial assumption that the 


follows: 


NEO is in an elliptical heliocentric orbit. Expressed mathematically, 


(zez (4.23) 


The NEO's true anomaly at impact (v NEO ) Jis fixed by the NEO's orbit. 
© / IMPA 


CT 


If this orbit crosses the Earth's orbit at perihelion, (re ) =0. If, on the other 
Lo ) mpacr 


hand, impact occurs at the NEO's aphelion, [v m =m radians. In between these 
76 IMPACT 


two extremes, the impact can happen at one of two separate true anomalies. Figure 4.4 
(above) demonstrates this fact. Thus, the true anomaly at impact in bounded as follows 


(expressed in radians): 


0c Us « 2x (4.24) 
ZS IMPACT 


To understand the bounds of the flight path angle ( y ), a graph of the relation 


between flight path angle and true anomaly is useful. This graph is presented in Figure 
4.5. 
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Flight Path Angle(gamma) vs. True Anomaly (nu) 
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Figure 4.5 Flight Path Angle of the NEO at impact relative to True Anomaly for various 
eccentricities 


One can see from the figure that the flight path angle varies between ln and 5 


Died. * (4.25) 


In conclusion, a study of the bounds of each parameter does not readily reveal a 


relation for V,. Thus, a plot of Vo may be more helpful. 


2: Graphical Display 


Since V., is a function of only two variables, it is possible to generate a three-dimensional 
plot of V4 against e and a. A small MATLAB m-file was written to generate such a 


plot. The m-file is listed in Appendix A. The plot is displayed in Figure 4.6. 
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V_infinity vs. a,e 
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Figure 4.6 V„ as a function of NEO eccentricity and semimajor axis 


The flat portion of this graph represents NEO orbits that do not intersect the - 
Earth's orbit. Note that for a given NEO eccentricity there 1s a semimajor axis that 


generates the maximum V. This is demonstrated graphically by the plot in Figure 4.7. 
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Figure 4.7 Variation of V„ for a fixed NEO eccentricity 


A second effect apparent in Figure 4.7 is the increase in Va as the NEO's 


eccentricity increases. 


3. Closed Form Solution 


As a final attempt to understand the geometry that defines V4, the symbolic 
manipulator Maple VTM was used to generate a closed form solution for the maximum 
value for V,. The Maple VTM code is included in Appendix B. 

The first step in the effort to find a closed form solution is to combine all of 
Equations 4.18 into one equation. Because of the assumption that the NEO is in an 
elliptical orbit about the Sun, it is possible to reduce the Maple V™ solution to the 


following equation: 
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© 
| a Y ek | 


(4.26) 
Next, the square of the magnitude of the velocity was calculated by adding the 
square of both components of the vector in Equation 4.26. In order to simplify the 
expression for the location of maximum Va, the function was left as the square of the 
magnitude. By eliminating the square root, the expression is simplified, but the final 


result will be the same. After simplification, the magnitude of V, reduces to 


ES. PR a D 
"^et of + 6 Ue Fay a “er SU Y, Y, 22 «(e en, 2a) Hai A, di 


V = 
rar -2a) 
% \% 


(4.27) 
Since the plots in the last section show that there is a semimajor axis (a) for each 
eccentricity (e) which generates a maximum Va, it will be most instructive to find an 
expression for maximum V,, with a fixed eccentricity. Thus, Equation 4.27 is 
differentiated with respect to the semimajor axis, with the eccentricity held constant. 


This derivative reduces to 
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(4.28) 


Solving Equation 4.28 for the semimajor axis gives 





(4.29) 


This relatively simple equation, considering the complexity of the derivation, is an 
important result. The equation determines the semimajor axis that will generate the 
maximum V. for a given eccentricity. 

Mathematically, there is no guarantee that Equation 4.29 is a maximum. Since 
only one derivative has been performed, all that is known is that Equation 4.29 is an 
extremum. In the Maple V™ procedure, inputting Equation 4.29 into the second 
derivative produced a very complicated result. Analysis of this result did not confirm 
that Equation 4.29 was a maximum. However, an eccentricity of 0.9 gave a semimajor 
axis of 1.739, which matches the maximum of the plot in Figure 4.7. This confirms that 


Equation 4.29 matches the graphical result and yields a maximum.. 


G. MAPPING TO A GEOCENTRIC ORBIT 


As stated before, two parameters are required to define a hyperbolic orbit: V. and 
the semiminor axis (b). The last sections went into great detail about determining V... 
However, knowledge of the NEO's heliocentric orbit is insufficient to uniquely fix the 
geocentric semiminor axis. 

In order to fix the semiminor axis, one would require exact knowledge of the 
location of the NEO, with respect to the Earth, as it crosses the Earth's sphere of 


influence. All that is known about the impact problem is the NEO's heliocentric orbital 
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parameters, and the prediction of a future impact with Earth. Knowledge of the 
heliocentric orbit only fixes the NEO's V, as it crosses the sphere of influence. However, 
there is no single semiminor axis that will define an impact based on V». 

.. A parameter used in the patched conic approximation for interplanetary arrival 
targeting can be applied directly to the impact problem. This parameter will determine a 
minimum semimajor axis required to avoid impact. The definition of this parameter is 
first based on the B-plane. Shown in Figure 3.1, the B-plane is perpendicular to the 
asymptote of the approach hyperbola, and placed at an infinite distance from the Earth 
(Brown, 1992, p. 116). 

In this Figure the semiminor axis (b) is a radius, drawn in the B-plane, from the 
location that the NEO pierces the B-plane to a line perpendicular to the B-plane that 
intersects the center of mass of the Earth. The impact radius (bj) is the semiminor axis 
that defines a hyperbolic orbit which is tangent to the surface of the Earth. Whenever the 
NEO passes inside the impact radius, it will impact the Earth. Thus, the impact radius 
defines a minimum safe semiminor axis. The impact radius is a function of V„ and can be 


calculated from the general hyperbolic relation (Brown, 1992, p. 28) 





D zd 


2 | (4.30) 


For the impact radius specifically about the Earth, 


2 
b=R, E E +1 (4.31) 
Gi = 


Since the previous analysis determined V. as a function of the NEO's heliocentric - 
orbital parameters, it is now possible to define the Impact Radius as a function of the 


NEO's orbital parameters. A plot of this relation is included in Figure 4.8. 
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Figure 4.8 Impact Radius as a function of the NEO's eccentricity and semimajor axis 


Note that, for clarity, the axes have been rotated from those used in Figure 4.6. 

Figure 4.8 confirms an important conclusion that can be reached from study of 
Equation 4.31: the impact radius increases significantly as V decreases. This is an 
important result when extended to analysis of the NEO impact problem. Hence, the 
greatest error in previous analysis (Knudson, 1995; Park, 1997; Elder, 1997; Ahrens 
1994; Meissinger, 1995; Solem, 1994) which neglected third-body effects arises when 
the NEO's orbit defines a small V4. These orbits are nearly circular, and close to the 


Earth (that is, with a — 1 AU). 
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V. APPLICATION OF THE PATCHED CONIC APPROXIMATION 


The preceding chapter established a theoretical foundation for applying a patched 
conic approximation to a three-body analysis of NEO hazard mitigation. Using this 
approach, a two-body analysis can be modified to include third-body effects. 

In their paper entitled "Minimum Delta-V For Deflecting Earth-Crossing 
Asteroids," Park, Elder, and Ross (Park, 1997) presented a trajectory optimization 
problem that was solved numerically using a sequential quadratic programming (SQP) 
method. The goal of the code was to determine an optimal AV to deflect an asteroid, but 
it did so using only a two-body approximation. The patched conic approximation has 


since been added to the code. The results are presented here. 


A. ORIGINAL TRAJECTORY OPTIMIZATION PROBLEM 


The performance index chosen to minimize AV is 


J = JAV? + AV? (5.1) 
where AV} is the component of AV which is parallel to the motion of the object and AV, 


is the component perpendicular to the motion. The constraints consist of the two-body 
equations that govern the motion of the NEO and the Earth about the Sun. Three terminal 
boundary conditions are set by the desired minimum miss distance. First, the miss 
distance at the final time must be equal to the minimum miss distance. Expressed 
mathematically 

EE (5.2) 
where R is the distance between the Earth and the NEO and R,riticat is the proposed miss 
distance from the Earth. The two additional boundary conditions come from the 


requirement that R be at a minimum when it reaches Roritical. Because R is continous and 


differentiable, this can be expressed as 
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R=0 (5.3) 
R>0 (5.4) 


B. THIRD BODY MODIFICATIONS TO THE ORIGINAL PROBLEM 


Third body effects can be added to the original formulation by modifying the 
boundary conditions in one of two ways. Both modifications are new definitions for 
Roritica as used in Equation 5.2. 

In the first modification, Reritica is redefined to be equal to the impact radius (in 
the B-plane) of the original NEO orbit. Once the original orbital parameters of the 
impacting NEO are determined, the user can input the desired miss distance and the 
modified code computes the impact radius as a function of the original NEO parameters. 


This can be expressed mathematically as 

R= M A original» SEPAral ion) =0 (5.5) 
where "separation" is the desired minimum separation between the NEO and the center of 
mass of the Earth. 

The second modification defines Reriticas as the impact radius of the NEO orbit 
after perturbation. Instead of defining R,itica; based on the original NEO orbital 
parameters, this modification recalculates the impact radius after the perturbation is 
applied and confirms that R,,itica 1S still equal to the impact radius. The second 


modification therefore changes the boundary condition to 


R-b(a,e,separation) = 0 (5.6) 
The two modifications produced slightly different results. The first modification 
changed only one calculation of b; and thus had little effect on the run time for the code. 
The second modification required several calculations of b; during execution of the code 
and thus increased the run time significantly. However, Equation 5.6 is a more direct 


statement of the required miss distance. This equation accounts for the fact that the 
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impulse which diverts the NEO will generate a new orbit. Since the impact radius is a 
function of the NEO's orbit, the perturbation will also generate a new impact radius. 
Equation 5.6 thus ensures that even the perturbed NEO will not pass within the impact 
radius. All results discussed below refer to the third-body modification expressed in 


Equation 5.6. 


C. SQP NON-LINEAR OPTIMIZATION ALGORITHM 


The trajectory optimization code written by Dr. Park is called Earth Defense 


against Asteroid Impact (EDAT). 


From the user's viewpoint, the code flows according to the diagram in Figure 5.1. 












Input: 
True Anomaly at impact 
2180 degrees? (y/n) 








Output to Plot: 
Orbits and 
Impact Point 













Input: 
Minimum Separation 
Distance (Earth Radii) 







Calculate: 
b,(a,e Minimum Separation) 













Input: 
Maximum impulse 
time (Asteroid Periods) 
AND 
Number of Steps 
per Period 









Output to Data File: 
Impulse Time, AV,, AV, 
t, Minimum AV 


Figure 5.1 EDAI Flow Diagram 
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At this early stage in the development of EDAI, the user 1s required to be very 
knowledgeable about the code's design. There are several parameters that must be 
adjusted, within the code, before each run. With some adjustments, discussed in the 
concluding chapter, the code can be much more user-friendly. 

As with all numerical optimization programs, the EDAI code is somewhat 
unstable with bad initial guesses. A critical parameter in the code is the initial guess of the 
minimum AV. If this initial guess is significantly different from the true value, the 
numerical optimization can become unstable, as is typical with numerical codes. 
Currently, there is no way to change this initial guess without getting into the code. 

Thus, if the user wishes to initialize anew EDAI run with NEO orbital parameters that 
are different from the previous run, he must input an accurate initial guess into the code. 
The initial guesses used for this analysis are included in Appendix C. 

Another instability comes from round-off error when trying to achieve a minimum 
separation of one Earth radius. Since the code uses AU as the distance unit, a minimum 
separation of 1 R, is extremely small. This is easily solved, due to the linearity of the 
problem, by running the code with a larger minimum separation (e.g., 10 R4) and dividing 
the result to achieve the desired separation (e.g., divide by 10). 

The results of each run are saved in a file which is named within the EDAI code. 
This technique allows for many kinds of analysis of the results. However, in order to 


change the file name one must make a modification directly to the code. 


D. RESULTS OF THIRD BODY MODIFICATIONS 


The impact radius (bj) is highly dependent on the orbital parameters of the NEO, 
as displayed in Figure 4.9. Since the third body modifications to the EDAI code make use 
of the impact radius, the differences in results are also highly dependent on the NEO's 


orbital parameters. 
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First, consider NEOs with a nearly circular orbit (e close to zero). In this case, the 
impact radius is relatively large. Thus, one would expect the difference between the 
results of the 2-body model and the 3-body model to be relatively large. Figure 5.2 is a 
plot of the minimum AV required to generate a 1 R, miss, as calculated by both the 2- 
body and 3-body codes, for a nearly circular orbit. (Note that, for standardization, all 
results will be presented at 1 R, minimum separation.) The x-axis is the time that the 
impulse is applied, in units of NEO periods before impact with the Earth. The y-axis is 
the minimum AV required in units of cm/s. 

Park (Park, 1997, p. 6) presents a thorough explanation of the results of the EDAI 
code. In summary, there are two effects apparent in the results. The first order variation 
in minimum AV required is dependent on the change in orbital elements of the NEO after 
perturbation. The second order variation shows local minima at the perihelion of the 
NEO orbit. These extrema demonstrate the commonly known effect that a AV applied at 
the perihelion has a greater efficiency in changing an orbit than when applied at any other 


place. 
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Figure 5.2 Minimum AV for I R, Separation (e=0.1, a=1.05 AU) 


For the original NEO orbit used in the plot for Figure 5.2, the impact radius ıs 
4.192 R,. Figure 5.3 shows that when the results of the three-body analysis are divided 
by the results of the two-body analysis the answer oscillates about a factor of 
approximately 4.195. This oscillation could be identical to the variation in the impact 


radius for various impulse times. However, this hypothesis has not yet been verified. 


Verification is suggested in the recommendations for future work. 
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Figure 5.3 AV(3-body)/AV(2-body) (e=0.1, a=1.05 AU) 


Figure 5.4 plots the difference between the 3-body and the 2-body models. The 
critical information in this plot is the increasing effect of the Earth's gravity as the impulse 
time decreases. This demonstrates the need to include the Earth's gravity in planning for 
short warning deflection scenarios. By not including this effect, the required minimum AV 


can be in error by more than an order of magnitude. 
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Figure 5.4 Difference Between 3-Body and 2-Body models (e=0.1, a=1.05 AU) 


Referring once again to Figure 4.9, it is apparent that the impact radius for highly 
eccentric NEOs is nearly one. Thus, one would expect that adding third body effects to 
the analysis of such NEOs would have little effect. Figure 5.5 is a plot of the minimum 
AV required for a NEO with eccentricity increased to e=0.9 and the semi-major axis 
remaining at a=1.05 AU. While it is nearly impossible to tell, both the 2-body and 3- 
body results are included in the plot. This very close correlation is a result of the fact 


that the impact radius for this original NEO orbit is 1.0593 Rg. 
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Figure 5.5 Minimum AV for 1 R, Separation (e=0.9, a=1.05 AU) 


While the effects of the third-body modification are significantly smaller than with 
the nearly circular orbit, the trends are the same with both NEO orbits. Figure 5.6 shows 


that the error increases with decreasing impulse time, just as in the first orbit. 
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Figure 5.6 Difference Between 3-Body and 2-Body models (e=0.9, a=1.05 AU) 


A critical benefit from adding the third body effects to the code is that the short 
warning results are now much more accurate. The final two plots are details of the 
minimum AV results for the two orbits for the timeframe of less than one orbital period. 
Note the sharp variation in the plot in Figure 5.8. The variation occurs when the NEO is 
near perihelion. As in the subsequent passes through perihelion, the AV requirement to 
perturb the NEO into a safe orbit is minimum at perihelion. Figure 5.8 shows that this 


applies during short warning time scenarios as well. 
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Figure 5.8 Short Warning Minimum AV for 1 R, Separation (e=0.9, a=1.05 AU) 


49 





E 


VL MITIGATION CAPABILITY OF CURRENT IMPULSE TECHNOLOGIES 


With a high degree of confidence in the new three-body minimum AV results, it is 
now possible to evaluate the potential of current technologies for deflection of NEOs. 
While Chapter II presented several ideas with potential for NEO deflection, the three- 
body analysis of this thesis was based on an assumption of an instantaneous AV (that is, 
an impulse). An impulse is practically impossible, and the only three methods discussed 
which approximate an impulse are kinetic impactors, a stand-off nuclear blast, and a 
surface nuclear blast. 

There are many variables to consider when analyzing the capability of an impulse 
technology, including the physical properties of the NEO, the mass and closing velocity 
of a kinetic impactor, or the explosive yield and stand-off distance of a nuclear 
interceptor. Ahrens and Harris (Ahrens, 1994, p. 897) present a detailed analysis of the 
effects of these variables. This analysis was combined with the minimum AV results from 


the last chapter to evaluate current deflection capabilities. 


A. KINETIC IMPACT 


A kinetic impactor imparts a AV on a NEO by the exchange of momentum from 
the collision of the impactor's mass and by creating a crater. The ejecta from the crater 
accelerates away from the NEO, thus causing a change in momentum. 

In their analysis, Ahrens and Harris obtain an expression for the mass ratio 


between the kinetic impactor (Mj) and the NEO (Mwygo) as follows (Ahrens, 1994, p. 


906) 
M; ` A / 1014222) | (6.1) 
My Yi | pi; X Y | 


Av — change in NEO velocity 





where 


v, = impact speed 

p = NEO density 

p, = kinetic impactor density 
Y = NEO material strength 

The numerator of Equation 6.1 gives the AV resulting from an inelastic collision, 
while the denominator accounts for the added momentum change from the crater ejecta. 

It should be noted that Ahrens and Harris differentiate between two different 
regimes in their analysis. For smaller NEOs, the strength of the NEO dominates the 
cratering mechanism. This effect is referred to as the strength regime. The gravity regime 
applies to larger NEOs where gravitational effects dominate cratering (Ahrens, 1994, p. 


904). The results were not significantly different, so the strength regime is used for the 


current analysis. 


B. STAND-OFF NUCLEAR BLAST 


A stand-off nuclear blast can impart a AV on a NEO by heating its surface until it 
expands and ablates away from the NEO. The ablated material is mass accelerating away 
from the NEO, thus causing amomentum change. Ahrens and Harris (Ahrens, 1994, p. 


912) conclude that the effect of a stand-off nuclear blast can be expressed as 


nAW 


Av = 0.1 T? 





(6.2) 
with 

n — efficiency of neutron production 

A = Geometrical efficiency factor (accounts for stand-off distance) 

W = nuclear explosive yield (expressed in Kt) 

D = NEO diameter (in km) 
Equation 6.1 expresses the NEO size in terms of mass, while Equation 6.2 uses 


diameter. For consistency, and for more direct application in the MATLAB m-file, 


Equation 6.2 was modified to give 


5» 


0.1InAWpz 
Mu "en — (6.3) 


e: SURFACE NUCLEAR BLAST 


The mechanism for momentum change with a surface nuclear blast is nearly 
identical to that of the kinetic impactor, except that the increased energy of the nuclear 
explosive creates a significantly larger crater. Since there is little experimental data about 
the effects of nuclear cratering in the strength regime, Ahrens and Harris present only a 
rough but educated estimate. The effect of a surface nuclear blast is therefore 
approximated by (Ahrens, 1994, p. 913) 

W = 4e - 5AvM zo (6.4) 


where Av is expressed in cm/s and Mygo in kg. 


D. ASSUMPTIONS 


Several assumptions are required in the analysis of mitigation technology 
capabilities. The physical properties of the NEO required by Equations 6.1 and 6.3 are 
the yield strength (Y) and the density. The yield strength can vary from 10° dyne/cm? for 
soft rock or ice to 10? dyne/cm? for hard rock (Ahrens, 1994, p. 906). For this analysis, 
a value of 10° dyne/cm? was used. The density of asteroids ranges from 2x10? kg/m? to 
5x10? kg/m? with a mean density of 3x10? kg/m? (Elder, 1997, p. 10). Comets are 
somewhat less dense than asteroids, with estimates of 1000 to 2000 kg/m? (Elder, 1997, 
p. 10). A mean asteroid density of 3x10? kg/m? was used for this analysis. 

Assumptions are also required for the mitigation technology. The goal of this 
analysis is to determine the maximum NEO size that can be diverted using a single 
mission with today's technology. It is always possible to send multiple diversion 
missions against a NEO threat, but a single mission baseline gives the most fundamental 


result. Rustan (Rustan, 1994, p. 1070) presents a summary of launch vehicle 
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performance which includes the useful payload that can be place on an interplanetary 
trajectory. He indicates that the Russian Energia vehicle can place a payload mass of 
18,040 kg into an interplanetary (or NEO intercepting) trajectory. Equation 6.1 implies 
that an impactor with high density would be most effective in diverting a NEO. Lead has 
high density and is inexpensive enough to be a reasonable candidate for use as a kinetic 
impactor. While the impactor would involve both guidance systems and concentrated 
mass, a fair approximation is to use the density of lead, which is 11 g/cm?. The closing 
velocity of the kinetic impactor is a function of the trajectory chosen to achieve intercept. 
For this analysis, a median closing velocity of 20 km/s was chosen. 

A search of unclassified documentation determined that the largest single nuclear 
warhead yield in a currently deployed strategic system is 24 Mt (Janes, 1996). This is 
the warhead used in the Russian SS-18 Mod 1 missile, appropriately nicknamed "Satan." 
Since no information about the warhead mass was available, it was assumed to be less 
than the 18,040 kg useful payload of the Energia launch vehicle. The neutron production 
efficiency (n) can vary from 0.03 to 0.3 (Elder, 1997, p. 43). A median value of n= 0.15 
was chosen for this analysis. The geometrical efficiency factor (A) is taken from a 
standoff distance of 0.4 NEO radii to give a value of A = 0.3 (Ahrens, 1994, p. 912). 

Finally, the results assume that the AV is applied in an optimal direction. 

All assumptions used for analysis of the capability of current mitigation 


technologies are summarized as follows: 
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NEO Physical Properties 
Yield strength Y = 10° dyne/cm? 


Density 9 = 3x10° kg/m? 


Kinetic Impactor 
Mass M; = 18,040 kg 
Density 9, = 11 g/cm? 
Closing Velocity vi = 20 km/s 
Nuclear 
Yield W = 24 Mt 
Neutron Production n=0.15 
Geometrical Efficiency A= 0.3 
Optimal Impulse Direction 


Single Diversion Mission. 





Table 6.1 Assumptions Used in Mitigation Capability Analysis 


E. TOUTATIS 


The asteroid Toutatis is an Earth-Crossing Asteroid that passes within 10 lunar 
distances every four years (Park, 1997, p. 9). Its low inclination makes it suitable for 
analysis with the planar methods used in this thesis. The orbital parameters of Toutatis 


are as follows: 


a = 2.5154 AU 
e = 0.6361 
i = 0.47° 


To apply the analysis to Toutatis, it must be assumed that i=0°. This condition 
is necessary to create an impact with the Earth, as well as to conform with the 


assumptions made in the patched conic analysis. 
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Toutatis takes the shape of two attached spheres which combine to give a 
diameter of approximately 4.3 km (Elder, 1997, p. 45). 

Figure 6.1 indicates the capability of the three impulse technologies against the 
Toutatis-type asteroid. The Figure shows a plot of the NEO's diameter versus the 
impulse time in Earth years. It should be noted that the critical properties that define the 
amount of energy required for mitigation is the minimum AV and the NEO mass. The 
NEO's density was assumed above. This assumption allows calculation of NEO 
diameter, assuming a spherical NEO. Since the NEO diameter is somewhat easier to 


identify with than its mass, the diameter is used in these plots. 
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Figure 6.1 Mitigation Capability Against NEOs in Toutatis-type Orbit 


It is apparent from Figure 6.1 that Toutatis-type NEOs cannot be deflected with a 
single kinetic impulse for many years. The diameter of Toutatis (4.3 km) lies 


comfortably within the capability of either nuclear option. However, Figure 6.2 shows 
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that a short warning impact from Toutatis cannot be deflected with any single mission. 
This Figure shows the deflection capability within one period of Toutatis. From this 
Figure, it is apparent that a surface nuclear impulse can deflect a Toutatis-type asteroid 
when the impulse is applied as close as 2 Earth years before impact. If the impulse is 


applied any closer than that, multiple missions will be required. 
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Figure 6.2 Short Warning Mitigation Capability Against NEOs in a Toutatis-like Orbit 


F. NEREUS 


For comparison, a second asteroid in a more circular orbit that Toutatis was 
analyzed. The orbital parameters of the asteroid Nereus meet the criteria, and are listed 
below (NASA Ames): 

a= 1.4897 AU 
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e = 0.3602 
i= 1.419 
D = 0.8 km 
Although this asteroid is significantly smaller than Toutatis, the long range plot in 
Figure 6.3 indicates that the kinetic impactor is still insufficient as early as 18 years 
before impact. However, it can be deflected with a single mission as close as 0.2 Earth 


years before impact, as indicated in Figure 6.4. 
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Figure 6.3 Mitigation Capability Against NEOs in Nereus-type Orbits 
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Figure 6.4. Short Warning Mitigation Capability Against NEOs in a Nereus-like Orbit 
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VIL SUMMARY AND RECOMMENDATIONS FOR FURTHER WORK 


A. SUMMARY 


The gravitational effects of the Earth have the strongest influence on minimum AV 
calculations for deflecting NEOs in nearly circular heliocentric orbits. Incorporating the 
Earth's gravitational effects can increase the minimum AV required for deflection by 
several times the value calculated using a two-body approximation. By combining a 
detailed understanding of the AV requirements for NEO deflection with the analysis of 
energy coupling by Ahrens and Harris (Ahrens, 1994, p. 897), an estimate of the 
capability of current impulsive technologies can be made. It is important to remember 
that the analysis presented in Chapter VI is not precise. Many approximations and 
assumptions were combined to provide an educated estimate of a single mission 
capability. Significant amounts of testing, including an actual NEO deflection, will be 
required to build upon and EI this analysis. Also, there were several assumptions 
made in the patched-conic analysis for simplification. While the results are accurate, and 
give good insight into the dynamics of a deflection mission, it must be remembered that 
the numbers are presented for NEOs in elliptical orbits that are co-planar with the Earth's 
orbit. The need to extend the analysis into the third dimension is discussed as a 


recommendation for further study. 


B. RECOMMENDATIONS FOR FURTHER STUDY 


1. Three-Body Truth Model 


A significant amount of work has now been accomplished towards understanding 
the orbital mechanics that govern a NEO deflection mission. The best way to test this 


analysis, barring an actual NEO deflection test mission, would be to build a computer 
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truth model. This computer code could model the three-body dynamics between the Sun, 
the NEO, and the Earth. A minimum AV from the above analysis can be applied in this 
model to confirm that the NEO is diverted as intended. 

The equations of motion for this code are a modified version of the solution to the 
classic circular restricted three-body problem. In the classic problem, two massive 
primaries are established in a circular orbit about their center of mass (Weisel, 1997, p. 
278). The third body orbits according to the gravitational influence of the two primaries. 
In the modified circular restricted three-body problem, one of the primaries (the Sun) is 
significantly more massive than both of the remaining bodies. As such, the second 
primary (the Earth) is established in a circular orbit not around the center of mass of the 


primaries, but the center of mass of the Sun. The geometry is shown in Figure 7.1. 





Figure 7.1 Truth Model Geometry 


The ij-coordinate frame rotates about the Sun's center of mass at a rate œ equal to 


the rotation rate of the Earth about the Sun. The Earth lies along the i -axis at a fixed 


distance x, . 
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First, consider the equations of motion as governed by Newton's Laws of Motion 


and Gravity (Wiesel, 1997, p. 35) 


N Gm.m. r, —r. 
mË =$ — 2 (7.1) 


Fi; je 





Applying Equation 7.1 to the three body problem shown in Figure 7.1 gives 
. Gmem Gm,m 
m,r- si) E) et 2:3 (7.2) 
z 3 


The zero in the first term of Equation 7.2 comes from the fact that the Sun is 


located at the inertial reference point. Dividing through by m, and performing some 








vector algebra gives 
x Gm Gm 
rn -- XE — = B (7.3) 
2 13 


The choice of r, rather than r, for the equations of motion now becomes apparent. 
By using r, the mass of the NEO cancels out while the mass of the two primaries (the 
Sun and the Earth) remains. Thus Equation 7.3 explains the gravitational force of the Sun 
and Earth applied to a NEO per unit mass. 


Some manipulation will provide the acceleration in terms of the coordinate system 








components 
7 Gm “ ay Gm ^ ^ 
A im (x,i+y.j)- iH (xi *ysj (7.4) 
2 3 
m m ^ m m ^ 
= AS, + "| d, +—y, i (7.5) 
ry n P n 


A useful result would be to reduce Equation 7.5 to two variables (x; and y?). This 


can be accomplished by use of some vector relations. 


At st (7.6) 


r, -r,-r, (7.7) 


Thus 
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^ 
* 


xj y,j=(x, —x,)i+(y,-y,)j (7.8) 


In the i-component, X» Is a variable but x; is a known fixed distance of 1 AU. 
The j-component reduces to y? because y; is defined to be equal to zero. Placing these 


results into Equation 7.5 gives 


(7:9) 





The two radii can be expressed in terms of x», y,, and the known variable x, as 


follows 


r, = fx ty (7.10) 


n= (x,-1AU) +y? (7.11) 
Equation 7.9 is the equation of motion of the NEO in the three-body coordinate 
system. Now it is necessary to investigate the dynamics of the system. Beginning with 
the definition of the r, vector, the velocity can be derived as follows: 
r, 2xityj (7.12) 
Kit nit ty] (7.13) 
The Transportation Theorem (Greenwood, 1988, p. 49) must be applied because 


the coordinate system is rotating about the k -axis. This Theorem states that 


izokxi 2a (7.14) 
and 
j=akxj=-ol (7.15) 
Substitution of Equations 7.14 and 7.15 into Equation 7.13 yields 
E, = xi x (oj) ye» (7i) (7.16) 
- (x, — y, )i * (y, * ax, )j (Tall 


The velocity can now be differentiated to give the acceleration 


. 
^ 
> 


r, = Ee - oy, Ji (x, - oy, Ji +(¥, tox); + ax, )j (7.18) 
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7 (X, -2ay, - o x, Ji (y, - 20x, 7 y, )j (7.20) 
Equation 7.20 defines the dynamics of the three-body problem in the rotating 
coordinate system. To convert to a state-space form useful in building the truth model, 
the components of Equations 7.9 and 7.20 are equated. This process yields the following 


final relations 


Im m | 
E as + E ) 


2 3 


A : n "m 
Y, * 2o, - ey, 2 -Gy, E +2) (mo 


2 2 
r Gar Fr 


ry = f(x, -1AU) +y? 


h hj 





With Equations 7.21, the motion of the NEO has been described using two second 
order differential equations. Because the equations are non-linear, the best solution is to 
convert to a state-space form and resolve them numerically on a computer. The truth 
model can then confirm that r, equals the desired minimum as calculated in the EDAI 


code. 


2: Three-Dimensional Analysis 


The assumption of co-planar orbits greatly simplified the above analysis, but it 
also severely limited the practicable application of the analysis. Generalizing the theory 


to three dimensions is critical in application to the majority of hazardous NEOs. 


3. Hyperbolic Orbital Analysis 


There are two points where a detailed hyperbolic analysis can be added to the 


NEO mitigation analysis. First, many NEOs can be established in a hyperbolic orbit 
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about the Sun before impact with the Earth. These orbits were not included in this thesis. 
Second, the effectiveness of an impulse applied within the sphere of influence of the 
Earth was not investigated in this thesis. Both areas are important for a thorough 


understanding of the NEO impact problem. 


4. EDAI Code Development 


The EDAI code has great potential to be very useful in designing a mitigation plan. 
As currently written, the code generates a minimum AV required to generate a given miss 
distance. By adding the analysis included in Chapter VI, the code could produce a more 
useful result, namely the ability of a proposed mission to effectively mitigate a known 
threat. In a fully developed code, the user should be able to input NEO's orbital 
parameters (including the inclination with respect to the ecliptic), the NEO's physical 
properties, and the properties of proposed mitigation efforts. The code would then 
output a plot of the mitigation capability against the diameter or mass of the threat-NEO, 
as in Figure 6.1. The user could then determine if the mitigation effort is sufficient to 
handle potential errors in estimates of the NEO's physical properties. 

As the first user not involved in the writing of the code, it would be constructive 
to provide suggestions for improvement. With some improvements and further testing, 
the code could be widely used in researching the hazard mitigation problem. 

The importance of an accurate initial guess was discussed in Chapter V. As 
currently written, a user is required to directly modify the code to input a new initial 
guess. A possible improvement would be to allow the user to input an initial guess. 
However, this would require a fairly well-educated guess, which is difficult to make 
without a thorough understanding of the problem. Another option for improvement 
would be to build a reference file of results from previous successful runs. This reference 
file could be accessed once the user has input the orbital parameters of interest to get an 


accurate initial guess. 
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Chapter V also discusses an instability induced by round-off error. The solution 
to this problem is to run the code with a larger minimum separation (e.g., 10 R,) and 


divide the result to achieve the desired separation (e.g., divide by 10). A more robust 
code could test for potential round-off instabilities and eliminate this complication for the 
user. 

It would be helpful to have the ability to name the output file during initialization 
while the code is actually running. This would eliminate the need to alter the code to 
rename the output file. Also, adding the NEO's orbital parameters and the required 
separation distance to the output file could help to avoid confusion between different 
files. 

A minor bug in the code was discovered when analyzing NEO orbits with a semi- 
major axis equal to 1 AU. In this unique case, the period of the NEO is identical to the 
period of the Earth. Thus, an impact would occur after each NEO period. The EDAI 
code, however, only calculates the minimum AV required to divert the NEO for a specific 
impact time. The code does not recognize that impact will occur after each orbital period. 
Thus, the results indicate that a AV applied several periods before impact will have a 
greater effect than a AV applied within an orbital period of impact. In reality, the results 


should be the same for each period. 


5. Billiards Scenario 


There is a fourth potential impulse concept that was not investigated in this 
analysis. Melosh , Nemchinov, and Zetzer (Melosh, 1994, p. 1115) suggest diverting a 
smaller NEO into a large threatening NEO. It could be useful to investigate this proposal 


using techniques developed in this thesis. 
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6. Verify Results 


The hypothesis discussed in Chapter V that the three-body results differ from the 


two-body results by a factor exactly equal to the impact radius needs to be verified. 
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APPENDIX A. V,, PLOT M-FILE 


%  Vasterth3.m 

% This m-file plots the magnitude of the velocity of an asteroid, relative 

% to the Earth, at the impact point. The variables are a_ast_sun and e ast sun. 
€  Assumptions are: 

% 

% - The Earth is in a circular orbit 

clear 


% First set some constants 


mu sun - 1.32712438e11; €$ km^3/s^2 
mu erth - 3.986005e5; % |km^3/s^2 
r erth sun - 1.4959787e8; % km 
R_sun = 6.96e5; % km 
a_min = (r_erth_sun + R_sun)/2; % km 
a_max = 10*r_erth_sun; % km 
Re = 6378; % km 


n max - 50; 


€ The first FOR-loop varies e from 0 to 1 at a step of 0.1 
for n = 1:(n max - 1) 


% The next (imbedded) FOR-loop varies a from a min to a max 
for step - 1:(n.max - 1) 
e(n,step) - n/n max; 
a(n,step) - (step - 1)*(a max - a.min)/50 - a min; 


€ Check that the orbital parameters define an orbit that intersects the Earth's 
€$ orbit. 
if a(n,step)*(1-e(n,step)) » r erth. sun 


V ast erth(n,step) = 0; 
bi(n,step) = 9; 


elseif a(n,step)*(l+e(n,step)) < r_erth_sun 


V_ast_erth(n,step) = 0; 
bi(n,step) = 0; 


else 


€ If the asteroid’s orbit will intersect the Earth's orbit, compute the true anomaly 
€ of the impact and the flight path angle at impact. 

nu imp erth(n,step) - acos((a(n,step)*(1-e(n,step)^2))/(e(n,step)*r erth sun) - 1/ 
e(n,step)); 

gamma imp(n,step) - atan((e(n,step)*sin(nu imp erth(n,step)))/(1 + e(n,step)* 
cos(nu. imp. erth(n,step)))); 
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% Then compute the magnitude of the velocity of the asteroid with respect to Earth. 
Vx1 ast erth(n,step) - sqrt((2*mu sun)/r .erth sun - mu sun/a(n,step))* 

cos(gamma_imp(n,step)); 
Vx_ast_erth(n,step) 
Vy_ast_erth(n,step) 

sin(gamma_imp(n,step)); 
V_ast_erth(n,step) = sqrt(Vx_ast_erth(n,step)A2 + Vy_ast_erth(n,step)A2); 
bi(n,step) » sqrt((2*mu.erth)/(Re*V ast erth(n,step)^2) + 1); % Re 
end 


Vxl ast erth(n,step) - sqrt(mu sun/r. erth sun); 
sqrt((2*mu. sun)/r . erth sun - mu_sun/a(n,step))* 


H H 


[Vmax(n),i] = max(V_ast_erth(n,:)); 


aVmax(n) = a(n,i)/r_erth_sun; 
eVmax(n) = n/n_max; 
[bimax(n),j] = max(bi(n,:)); 
abimax(n) = a(n,j)/r_erth_sun; 
ebimax(n) = n/n_max; 
end 

end 

figure(1) 
colormap('white') 
aplot - a/r.erth sun; 
surf(aplot,e,V ast erth) 
grid on 
xlabel('a (AU)') 
ylabel('e') 
zlabel('V inf in ity (kmn/s)') 
title('V inf in i.t y vs. a,e') 

figure(2) 

l colormap('white') 

surfl(aplot,e,bi) 
grid on 
xlabel('a (AU)') 
ylabel('e') 


zlabel('b_i (Re)') 
title('Impact Radius (b. i) vs. a,e') 
view(-37.5-125, 30) 


answer - input('Generate 2-D impact radius plot? (y/n) ','s'); 
while answer=='y'; 

e interest - input('At which eccentricity? '); 

index = round(e interest*n. max); 

figure(3) 

plot(aplot(index,:),bi(index,:)) 

grid on 

xlabel ('a (AU)') 

ylabel (‘bi (Re)') 

e title-num2str(e.interest); 


title(['bi vs. a for e = ',e title]) 
print ans - input('Print 2-D plot? (y/n) ','s'); 
if print ans — 'y' 
figure(3) 
print 
end 
answer = input('Generate 2-D impact radius plot? (y/n) ','s'); 


end 
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APPENDIX B. MAPLE VTM ANALYSIS FOR V, wax 


File name: vbnd9.mws 
Last Modified: 5 August 1997 


Written by: LT Scott Porter 


Purpose: 1. Determine if the impact velocity between the Earth and an asteroid in an elliptical 


. Orbit has a maximum. 


2. If a maximum exists, find a function of an asteoid's semi-major axis and it's 
eccentricity which results in the maximum impact velocity between the Earth and an asteroid. 


emp o em u^ amb @@ amb ue Gub HD cu we Gub Gub em em om em em ar a uw gr em emt em UD UD AUD mb em emt em em em SF SS SFT SF SF SS SSF SF om em e em em SF cum Gub Gub AED: GUD emt em om em gz em om em om az e em em em om em ar om em om em om az er ape Em aue pe em rz cus dem S88 SF S| ER sum gr ar om emm emt em em. emt em em em em em em em em em em em 


This analysis uses the prinicples of orbital mechanics, vector algebra, and simple calculus to determine 
if the relative velocity between an asteroid and the Earth at impact has a maximum. The asteroid is 
assumed to be in an elliptical orbit which is co-planar with the Earth's orbit. The Earth is assumed to 


| be in a circular orbit. 


, 2 restart: with(linalg): readlib(isolate): 


Warning, new definition for norm 
Warning, new definition for trace 


Express the magnitude of the velocity of the Earth and the asteroid. 


> v[E] (a) :=sqrt (mu[sun]/r(E,sun]); 
> v[A, sun] (a) :=sqrt (2*mu[sun] /r[E,sun] - mu[sun]/a); 


a) = 





V4 f BEN a ) C 





Convert these velocities to vector form. The coordinate system is defined with the X-Y plane in the 
orbital plane, the origin at the center of mass of the Earth, and the Y-axis continuosly pointing to the 
center of mass ofthe Sun. Thus, the velocity of the Earth is always fixed along the X-axis. The 
velocity of the asteroid is defined by the flight path angle (Gamma(a)), the angle between a tangent to a 
circle at the asteroid's radius and the velocity vector. 


> VIE, sunl (a) :=vector(2,[v[E] (a) ,0]) ; 


> V[A, sun] (a) :=vector(2,[v[A,sun] (a) *cos (Gamma (a) ) , v[A, sun] (a) *sin(G 
amma (a))]); 


Fl 





| 
| 


Hun 
Fe S SE S 0 


] E. sun 





d 
H cun U cun : | 
2 [ae = ae a 


Cl 









En A m 
mcos [gas 
a 


| Pan T ^ 


pS E . 
E. sun ën 


Using vector algebra. define the velocity of the asteroid with respect to the Earth (V[A,E]). V[A,E] is 
defined as the velocity of the asteroid with respect to the sun (V[A,sun]) minus the velocity of the 
Earth with respect to the sun (V[E,sun]). 


» V[A,E] (a) :zevalm(V[AÀ,sun] (a) -V[E,sun] (a) ): 


3 Heun Heun Usun H sun M oun MT 
SE SEHR DECUS 2 = snl CoB 
! E. sun a F E. sun / E. sun e 


Now determine the orbital position and flight path angle of the asteroid at Earth impact. 










lira. | 


> nu[I](a):-arccos((a*(1-e^2)/(e*r(E,sun])) - 1/e); 
> Gammal (a) :=arctan(e*sin(nu[I] (a)) / (1+e*cos (nu[I] (a)))); 


( alte) - 
WER E 


e E sun € 
e ] | 


[I l(a):sarctan| a aa 
(2d — e`) - 
lela oe 





> Gammai (a) :=simplify("): 





S 
EE 2 
EE +247.) ace 





> Pe sun 


e PE. sun 
I l(a) := -arct > 
a(-l +e’) 


Substitute the impact orbital parameters to determine the velocity of the asteroid with respect to the 
Earth at impact. 


» V[A,E] (a) :2subs (Gamma (a) Gammal (a) , V[A,E] (a)); 
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í 2 NES : 2 
EE E n LA E 0) 















| E > 2 a sun 
H sun H sun | e Fg . sun 
— - SE —arctan| 
VE. sun e \ al-l +e) 











/ > a a 2 E 
Geren (ea, +2ar- =a ) 
E = ^ 2 DE sún 
8 € FE. sun 
sin| —arct : 
a(-1+e”) 


Now go through several steps to simplify the expression. 


> ViA,E] (a) :=vector(2, [simplify (V[A,E] (a) [1]) ,samplify(V[A,E] (a) [2]) 
1); 












Hun i sun i a) e sun (Ze sun +2 a) 





(= etic!) a 





a Pe sun 





te sun +2 a) 






Pe sun 


(-1 rena 





> SE 2 > 
DE td) (el +e ) (a € ~ FE, sun Zar, a ) 








» } E. sun 


a Fg sun 


2 
€ Fg sun 





VE sun ice sun +2 a) 


(-1 rea an | 
> V2[A,E] (a) :=vector (2, [combine (V[A,E] (a) [1]) ‚combine (V[A,E] (a) [2])] 












| n 

| Asun rn, +2 a) V Ce FE. sun = 2 Heun a) (-a ae a e”) 

i ; ATE, sun a pa e 

| V2, ga) = Ba 


Er sun (Te. sun + 2 a) 


(-l+e)a 


E T > 4 4,2 ` 3 ae 
Hun a e FE. sun E Asun @ e PE sun Y Hus d € — Hun lE sun un PE. sun a 13 H sun PE sun @ 


JS 


D IE 
SE l 4 » : a 
+12 Hm Tre —-du ada) AE 


sun ‚sun 


SENS [A, E] (a) :=vector (2, [expand (V2[A,E] (a) [1]) ‚expand (V2[A,E] (a) [2])] 


Eee —--- 


"9 





eh a 





> SR ^ à 3 
= : > . LET D Ce 2 
H cun D sun ar H iun N sun deti H sun a = H sun a e 


4 








de F d r 
E. sun E. sun E. sun E. sun > 3 
3 ms > ? 2 7 (-a «a e) 
(-1+e)a (-1+e)a (-l-e)a (-1+e)a 


3 


A 4 2 s 2 v 6 2 
MT d e PE. sun E sun ae a sun d Asun ae ma Houn Fe sun Së H sun PE. sun a 
= 3 4 2 Ur? 
S - st ^ 
—-13n feu. 0 tl2gg no D / (a( EE M 


> V4(A,E] (a) :=vector(2, [simplify (V3[(A,E] (a) (1]) ,simplify (V3[A,E] (a) [ 
2a, 





H sun br sun +2 a) H sun E. sun t2 a) 


» 


a + 








PA 


a Pe sun 


a FE, sun 


A (-1 +e") a’ are +2 2| 


2.02 2 2 2 
Aa e In Zar on mi e 





a (Stee de a) Zr sun 


Find the square ofthe magnitude ofthe velocity. Using the square ofthe magnitude, rather than just 
the magnitude, simplifies the algebra but gives the same final result. The maximum of a function occurs 
at the same place as the maximum of it's square. 
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— 


EN ua3g Au -USDIA Ei(atll)^2 - V3I[A,E] (a) [21^2: 





S 
—H un Ge sun aoc H cun E um d cs dis Houn E E Houn C 









> 


bo sun WC sun NET sun Bn sun ^ a 9 
RE  — = BD lee) 
(-l+e)a (-l+e)a (-1+e)a (-1+e)a 

E 
SNE 2 E > 4 Bras 
e Henn a e ME sun 4 L un a e e .sun ze Houn f em a Hyun lg aun O Hyun SE ae 


3 


3 SÉ 
ex a 3 2 . = E 
= 13 Asun Pe. sun a + (2 Hoon eu FE sun 4 Houn a ) J (a (~ E. sun +2 a) ae sun ) 
Simplify the expression. 


> Vmag[A,E] (a) :=simplify("): 


3 


r ET = 2 ’ 
Vmag, (a) 17) FE, sun H sun ar 6 H sun a + Hun lE sun 








Houn (= we >? a) 


+2 


(ON "Hu el + e ) a En sun Y E a) A Lag a (See +2 a)) 





a e sun 


> combine("); 


2 - > 
=D FR sun H sun a+ 6 Usun a + H sun B sun *2a( 
3 4 ? 3 2 3 4 7? 


= = d wes DE = CEU ne SX 
(Houn a FE, sun E44 a Bud Ho a — Yun de Bm all, ae m Et a 


l 4 


SEH en) ee) 
factor ("); 


= 
A 


2 
= d 2 r 
> Ca 6 ees yh T H an TE en 


*t2a4- sun esos O 2209. PAD a ( le sun + 2.4)) 


> collect(", S 


(6 Moun PE. sun a E (—5 H sun D +2 V SH a (e EE ) (e +1) (SA +2 a) PE. sun ) a 


3 2 
+ Asun FE. sun ) e Gr a S sn +2a)) 
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| 
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— 
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Es 


| 
| 


i 


J 
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> Vmagsara[A,El (a) :=sortí("): 


2 2 H : ? 
| magsqáqi d, al A= 6 cl H cun SR sun "1 H cun E sun 


+(-3u 


i Sun 





? 3 
B — — 2 2 zes. e: ue en E 
UE. sun +2 (€ |) (e + l ) ( ie uU UO Ma TE. sun ya) 7 (( 2 E a 


> 


be sun ) 


Differentiate the function for the magnitude (squared) and set the result equal to zero to find the 
maximum. 
> Vmax1l[A,E] (a) :=simplify (diff (Vmagqsqrd[A,E] (a) ,a)=0) ; 


E f A ED 3 ^ SN > . _ . dM 
max] El a) : ( = H sun are Ll de +a H oun om Hun ar sun de 


sun 





4 > 
d = > 
= — 3 — $ . e a 
(e l Ie era Vp sane OH ërem len / (een al 


2 ES sun 





Ia l) (e zl ( 2d— Pe sun ) a Bon PE. sun al 
> Vtest[A,E] (a) :=simplifv(dıff(Vmaxl1[A,E] (a) ‚,a)): 


4 


F vum ie E 3 E D 2 = a e 2 g 
Vtest , (a) B ^ ( 2 Heun qaot H sun aeta H sun VE. sun H sun a sun ae 


leer -(e-1)(e*1)(2a-rg,u) aM. ream) Bas Io 





-(e = ) (e SR ) (2 a — FE. sun ) a Ham FE. sun Geer -0 


Solve this function for the semi-major axis. This result will determine, for a given eccentricity (e), the 
semi-major axis which will generate the maximum relative velocity between the asteroid and the Earth 
at impact. 


> isolate (Vmax1[A,E] (a) ,a); 


Rca Vee 
(Sl tem) ) PE. sun 


a 3 
=1+e” 


> subs (a=amax,"); 


320813 
See) 
amax = 


VE. sun 
2 
=l+e 
> assign("); 


To further simplify the relation, use r[E,sun]-1 AU. 


> aAU:=subs (r[E,sun]=1,amax) ; 
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M E 


(-(-1 me) ) 
BU en 

` -] +e 

ME SS TRATA EDS -Isubs(acaAU vVtest[(A E] (a)); 


( 


y y 


2 ^ = = 
A Bun Houn € ( E = mee ) ) H An e sun 
) 01) ———— EE 


| 

Hesi2 (a) i= — > Ee ae : 

pa — — e 2 
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wo. 
M oun D sun A Ser ED ) e 
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(-(-l*-e€)) , 23 10/3 2 
Waller) > e T FE, sun lern) Mun PE. sun 
2: e 3 
xu 


(-1+e%) Mi 





2 






meee TS 
-l-e)) AE > 
(em Die) 2 > — FE, sun (El +e) ) Asun FE. sun 
-l +e 
m > FE, sun 
-l+e 


Sch) 
Szsimplify("); 


l sr. vol e +2 Asun e -2 HE TE poles el +e) Po H sun Ken) 
2 FE, sun y 7o] 


Ee 2 Ben, S 
ee +e) ) es: e^) (Cl Te) ) Usun FE. sun 
| %l DI Tg ILIO ‘G 

ee, 
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> subs (r[E,sunj=1,"); 


moo a nr 
l Hm (4 % | 9 ONT cud vol Bcc +e) ) Hes) 


i 
2 "PE 


goles EE 2221/3 
(2(4-1+e)) -I+re) (-A-I+e)) Hyun 
ee nn 


t3 


e EE 


| Hu TN Dä vg 1 3 Hun  — 3 Han EVAL LHe) Wan) 
2 A 961 

TRIS Tes Ea O D 

(AMELIE A 


Yl i= - A 
i (see) 


> eval (subs (mujsunj=1,"}); 
| 44/94] e* + e722 = 4/9 ERES 
E 4| 961 


s SE N, 
(-2(-1-1+e)) -1+e) (-(-l+e) ) 


= 0 


%l :=- 5 
(ee 


= Smpiscr um 
KENE ome my 0 
2 4 961 
el EE „2 “a 
(-2(-(-1+e)) —-I+e) (-(-l+e) ) 
%1 == 7 er 
(-l +e) 


i End of analysis. 
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APPENDIX C. INITIAL GUESSES FOR EDAI CODE 


As with any numerical optimization algorithm, the EDAI code requires an 
accurate initial guess of the results. An inaccurate guess will result in an instability. For 
this analysis, the code was first run with a known solution used to provide an initial 
guess. Then the NEO's orbital parameters were altered slightly, such that the initial guess 


from the previous run was still close enough to the correct answer to provide a stable run. 





This process was continued until a database of initial guesses was built. The initial 


guesses are provided here. 


Canonical units are used for the AV and ts guesses. The units are defined as 


follows: 


DU = Distance Unit = Astronomical Unit (AU) = 1.495978e8 km 


; l 
TU = Time Unıt = — year 
2x 


Minimum 
a(AU) i Separation (R,) 


1.062e-7 
4.963e-5 
-1.581e-5 
-8.548e-5 
-1.4417e-4 

-7.926e-4 
-6.057e-4 
-4.221e-4 
-5.087e-4 
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AV, (DU/TU) | AV , (DU/TU) : t: (TU) 


4.142e-4 
4.157e-4 
3.999e-4 
3.575e-4 
2.769e-4 
8.303e-4 
4.340e-4 
3.440e-4 
7.218e-4 


2.843e-4 
1.336e-4 
8.856e-5 


! 6910e-5 


5.543e-5 
1.674e-4 
3.446e-7 


| -133le-4 





| -4.042e-4 
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